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Abstract. We analyze the emergence of Shapiro resonances in tunnel- coupled Bose- 
Einstein condensates, realizing a bosonic Josephson junction. Our analysis is based 
on an experimentally relevant implementation using magnetic double well potentials 
on an atomchip. In this configuration the potential bias (implementing the junction 
voltage) and the potential barrier (realizing the Josephson link) are intrinsically 
coupled. We show that the dynamically driven system exhibits significantly enhanced 
Shapiro resonances which will facilitate experimental observation. To describe the 
systems response to the dynamic drive we compare a single-mode Gross-Pitaevskii 
(GP) description, an improved two-mode (TM) model and the self-consistent multi- 
configurational time dependent Hartree for Bosons (MCTDHB) method. We show that 
in the case of significant atom-atom interactions or strong driving, the spatial dynamics 
of the involved modes have to be taken into account, and only the MCTDHB method 
allows reliable predictions. 
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1. Introduction 

Coherent tunnelling dynamics of macroscopic many-body quantum states through a 
classically forbidden barrier is one of the most striking manifestations of quantum 
physics. In the presence of an external drive, such systems exhibit the Shapiro- effect^ 
also known as photon- assisted tunneling [K |2]. Shapiro resonances have been studied 
extensively in the context of solid-state superconducting Josephson junctions where a 
resonant modulation of the energy bias leads to a DC current across the junction [3j. 
Such resonances can be used to exactly quantify an applied DC potential difference 
across the junction and are nowadays used to implement a voltage standard. 

A similar effect occurs in tunnel-coupled Bose-Einstein condensates (BEC) in 
double well potentials, where a resonant modulation of the bias energy lets the 
condensates oscillate between both wells [2j. Such dynamics could be used for an 
accurate measurement of the chemical potential difference across the tunnel barrier [4J . 
Despite the close analogy to the superconducting system, the dynamics of the hosonic 
Josephson junction are strongly influenced by atom-atom interactions [5j . Additionally 
the presence of a finite atom reservoir inhibits the observation of a true DC atom flux. 

In recent years much effort has been devoted to the theoretical and experimental 
study of the tunnelling dynamics of driven Bose-Einstein condensates [6j. In optical 
lattice experiments, arrays of bosonic Josephson junctions have been realized [7]. 
Shaking of the lattice allowed to control the effective tunnel coupling [8j and the 
superfluid to Mott-insulator transition [9j. This effect was also used for dynamical 
localization [10]. The direct analogy to the superconducting Josephson junction has 
been realized in an all-optical double well potential. Here the influence of atom- 
atom interactions has been evidenced, leading to new dynamical modes [11] and sub- 
poissonian number statistics [12j. 

In this work we analyze the Shapiro effect for magnetic double well potentials 
on atomchips in view of a recent experimental realization of a bosonic Josephson 
junction We focus on the experimentally relevant situation where a modulation 

of the energy bias between the left and right BEC is accompanied by a concurrent 
modulation of the tunnel coupling. As we show, this leads to a significant enhancement 
of the Shapiro resonances compared to the conventional driving [2] , which will facilitate 
experimental observation. 

Shapiro resonances have been theoretically analyzed in [lH [2] within a two-mode 
(TM) model employing static condensate wave functions [T5l[l6]. For the configuration 
investigated in this work it is important to take into account the spatial dynamics, 
as transverse excitations might occur due to the external drive, which have a strong 
influence on the tunneling dynamics. We therefore employ the multi-configurational 
time dependent Hartree equations for Bosons (MCTDHB) method [El [IE], which is a 
self-consistent, essentially exact framework, taking into account the full dynamics of the 
system. In previous work, the MCTDHB method applied to bosonic Josephson junctions 
allowed to discover dynamics very different to that predicted by the Gross-Pitaevskii 
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equation and the two-mode model, such as an inverse sphtting regime [19], the decay of 
macroscopic quantum self-trapping [20j , and violation of time-reversal symmetry when 
switching interactions from repulsive to attractive in bosonic Josephson junctions [21j. 
In this work we study another case, where for typical interaction strengths the Shapiro 
resonances predicted from MCTDHB appreciably differ from the simpler models. 

The paper is organized as follows: We start with the physical system and its 
dynamic description in section[2j In section [2?T] we introduce the realization of the double 
well potential and how we implement the driving, and in section |2.2| the three models 
which will be used and compared in this work: the Gross-Pitaevski (GP) equation, 
the two-mode model, and MCTDHB. In section [3] the non-interacting driven system is 
discussed, and results which demonstrate enhanced Shapiro resonances are shown. It 
includes also a short discussion of how to choose experimental parameters on order to 
observe the effect. Results for the interacting case are presented in section [4j First we 
show how the resonances shift due to interactions in section 14. 1[ A discussion of the 
different dynamical models, as well as the damping of the number imbalance is finally 
discussed in section 14. 2[ 



2. Physical system 

2.1. Double well potential 

The following analysis is motivated by a recent experimental realization of a bosonic 
Josephson junction on an atomchip p^. A symmetric double well potential is generated 
by a combination of static and oscillating magnetic near fields, making use of dressed 
adiabatic states [22l [23l [24j [|| The system consists of two elongated traps with strong 
atomic confinement {uj^ ^ 27r • 2 kHz) in the x and y (transverse) directions, and a very 
weak confinement along the z (longitudinal) direction with (jJx,z/^z ^100. Tunnelling 
dynamics takes place along the coupling direction x which connects the two potential 
minima of the double well. The system is assumed to remain in the (many-particle) 
ground state in the two orthogonal directions y and z which do not contribute to the 
dynamics. The distance d between the double well minima and the height of the tunnel 
barrier are adjusted by controlling the amplitude of the oscillating field component which 
is in the radio frequency (RF) range. This RF amplitude represents the dynamical 
control parameter \{t) used to drive the Shapiro resonances. 

To implement a potential difference ("voltage") between both wells, the double well 
can be tilted in space (making use of gravity) [221, [25] . Further possibilities are a spatial 
inhomogeneity in the amplitude of the RF field or employing additional forces like DC 
electric near fields [26j or local optical dipole potentials [lU [27] . The resulting potential 
is given by equation 10 of [23] and can be approximately described by a symmetric 
fourth order polynom V\{x) for the double well with an additional linear gradient g that 

X The considerations made in this work are equally valid for double well systems based on solely static 
magnetic fields or on optical dipole potentials. 
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Figure 1. Q and AE versus A, calculated using eigenstates of the Schrodinger 
equation. Upper right inset: Driving geometry. Lower left inset: When the splitting 
parameter A (black lines) is driven periodically, both Q (violet lines) and (green 
lines) oscillate. Solid line: first resonance n = 1 for the non-interacting system, with 
parameters as in text. Dashed line: second resonance n = 2 (with larger amplitude). 

implements the tilt: 

Vx^g{x) = Vx{x)+g-x 

= ci(A) [x - d{X)/2f + C2(A) [x - d{X)/2]^ + g-x. (1) 

Note that the control parameter A (RF amplitude) affects both the separation of the 
double well minima (rf(A)) and the steepness of the confinement (expressed by the 
coefficients Cid{X) and C2rf(A)). By construction, a larger separation coincides with a 
higher potential barrier, both reducing coupling, making A the most sensitive parameter 
in the system. 

The most straightforward implementation of a Shapiro experiment would be 
to dynamically modulate the potential gradient g. However in a purely magnetic 
implementation of the double well as in [13j this would also result in a significant spatial 
displacement of the trap minima, which might lead to excitations and uncontrolled 
dynamics. We therefore focus the analysis on the effects of a periodic modulation of 
the RF amplitude (splitting parameter A(t)) keeping g fixed. This results in double well 
minima moving along a potential slope as indicated in figure [l] (upper right inset). As 
we will show, modulating both the tunnel coupling and the potential difference leads to 
a significant enhancement of the observed Shapiro resonances. 

2.2. Model 

2.2.1. Many-body Hamiltonian In this work we consider the dynamics in the splitting 
direction x, and neglect any dynamics in the other directions y and z [28j. The dynamics 
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of the atoms in the double weU is then governed by the many-body Hamiltonian [291 [30 



¥{x)h^{x)i!{x) + ^i}\x)¥{x)i}{x)i!{x) 



(2) 



where ^(x) is the bosonic field operator, and h^{x) = — + Vx^g{x) is the bare 
Hamiltonian. Uq is the eflFective ID interaction strength, obtained by integrating out the 
other spatial directions for the ground state j28j. The ID approximation is reasonable 
since typically the dynamics in the other directions decouples to a good extent whenever 
the 3D-Potential factorizes, i.e., V{r) = Vx{x) + Vy{y) + Vz{z). 

Direct solution of the dynamics due to this Hamiltonian is not possible for more 
than a few atoms, and thus we have to use approximation schemes. In this paper we 
will employ and compare several types of approximations which we will outline in the 
following. They can be considered as diflFerent ways of restricting the field operator ^^(x) 
to a small number of modes. 



2.2.2. Gross-Pitaevskii (GP) equation The mean field dynamics [29] is obtained by 
restricting the field operator to a single mode function, ^[x] = ao(/>o(^). The Heisenberg 
equation of motion for ^(x), obtained from equation [2| yields the GP equation 

Z(/)o(x, t) = (h\x) + Uo{N- 1) |(/)o(x, t) (/)o(x, t) . (3) 

For a BEG in a double well, the GP wave function can be written as 00 (x) = 
[(/)l(x) + (()R{x)]/^/2^ where {(I)r{x)) is localized in the left (right) well, and thus 

the many-body wave function reads 




The GP equation assumes that the number distribution between the left and right 
localized condensate has binomial number fiuctuations ^/N (perfectly "coherent 
state" [31j). Since for the initial state we consider in this work always a condensate 
localized in one well, the GP equation will well describe the short-time-dynamics. At 
longer times, however, a single-mode description will not be valid anymore. 

2.2.3. Two-mode (TM) model In a double well it is more natural to use a basis which 
comprises two modes instead of one: ^^^(x) = aL(j)^{x) + aR(j))^{x). Here, the localized 
mode functions (j)^{x) and 0^(x) are depend on the control parameter A. oll (a^) 
annihilates (creates) an atom in the left well, and similar for clr (a^). 

There exist several schemes for the choice of the modes (j)\{x) and 0^(x). The 
simplest one uses superpositions of the two lowest lying eigenstates (/>^(x) and (/>e(x) of 
the single-particle Schrodinger equation of the symmetric potential Va(x), which have 
gerade and ungerade symmetry: (p^^ix) = [(j)g{x)±(j)^{x)]/^/2 The Hamiltonian 

reads in terms of pseudo-spin operators: 

H = -n{t)X + ^E{t)J, + 2i^{t)Jl . (5) 
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Hereby, Jz = |(S^a^ — a^a^) measures the atom number difference between left and 
right weU, and Jx = K^l^i^ + ^i^^L) proniotes an atom from the left to the right well and 
vice versa. Equation [5] depends on the generic parameters AE{t) and /^(t), which 
denote the tunnel coupling, energy bias, and nonlinear interaction energy, respectively. 
They are given as 

n{t) = - J dxf£^{x)h\x)(t)^^{x) + h.c. , (6) 
AE{t) = J dx(l)"£^{x)h^{x)(^^{x) - J dxfjl^{x)h^{x)(^^^{x) , 

^0 f J^IaX /^m4 



The time-dependence in these parameters is due to a time-dependent control A(t). 
Tunnel coupling and energy bias between left and right condensates are shown in figure [T] 
over a wide range of A from the unsplit to the split case. [§] While the tunnel coupling 
decreases with increased splitting, the energy bias becomes larger. Typically, we choose 
a constant ^ t/o/2, which is a good approximation because this value depends less 
crucial on the shape of the modes as compared to Q{t) and AE{t). 

An improved two-mode model [32] can be obtained by using the first and second self- 
consistent states of the GP equation and (/)f^'^{x) with self-consistent energies 
E^^'^ and respectively. The tunnel coupling of the improved model f]^^^ is given 
by the energy diflFerence E^^'^ — E^^''^ and a shift due to interactions. The parameters 
for the improved model are then given by 

^^'\t) = ^GP-A _ ^GP,x _ u,N/2(^J dx\(t>':''''\^ - J dx\(t>f^''\^^ , (7) 

AE('\t) = j dx[(l>'l^'Yh'<t>T'' - j dx[ff^Yh^<f/^^ . (8) 
In the following we always use the improved two-mode model, unless stated otherwise. 



2.2.4' Two-mode model with self- consistent orbitals and occupations A dynamical 
description including also the spatial dynamics can be obtained by the multi- 
configurational time dependent Hartree equations for Bosons (MCTDHB) [17j, which 
represents a framework where the two modes are included self-consistently. Using two 
time-dependent orbitals has the crucial advantage that transverse spatial excitations 
can be included, in contrast to a model with two fixed orbitals. 

A full discussion of the working equations of MCTDHB [TTl [28j can be found 
elsewhere. We only sketch here the main ideas of the method. The ansatz for time- 
dependent modes reads ^^^(x) = aL{t)(j)\{x^ t) + aR{t)(j)\{x^ t) for the case of two modes, 

§ We use units where = 1, mass of a ^^Rb atom m = 1, and atom energy and time is scaled by the 
confinement length a^o — \J^I {rn^ho) — 1 l^rn and energy hjJi^o of a harmonic oscillator. It follows 
that the units of time and energy are, respectively, 1/hwho = 1.37 ms and hw^o = 27r • 116.26 Hz. 
II Other, more subtle corrections to the standard two- mode model (e.g., to /i:), as suggested by Bergeman 
et at. [32], will not be used here since it does not lead to significant improvements. 
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although in principle the number of orbitals can be chosen at willjj] MCTDHB then 
provides a way to determine the "best possible" shape of the orbitals at a given time. 
This is achieved by formulating an action integral based on the above ansatz, and then 
using a variational principle. For the orbitals one obtains nonlinear equations, and for 
the number part two-mode equations similar to the TM model. For the orbitals one 
obtains nonlinear equations, and for the number part two-mode equations similar to the 
TM model. Most importantly, the orbital part and the number part are coupled and 
thus have to be solved self-consistently. 

An important quantity is the one-body reduced density (^^'''(x)^^(x)) [35]. It can 
be diagonalized in terms of natural orbitals and natural occupations (i = 1, 2): 

2 

p{x) = Y,pmx)\' ■ (9) 

i=l 

Whenever one natural occupation dominates, we have a BEC [36] and one orbital 
suffices. MCTDHB then coincides with the GP equation. Whenever several eigenstates 
are finite, we have an m-fold fragmented BEC, and m-orbitals should be used. Whenever 
MCTDHB-calculations converge when increasing the number of orbitals, the results can 
be considered as an exact solution of the many-body Schrodinger equation. 

In this work we measure the degree of fragmentation by the difference of the 
population of the natural orbitals. (pn — P22)/N = ±1 corresponds to a single BEC, 
whereas (pn — P22)/N = corresponds to a fully fragmented system. 

2.3. Observable 

The Shapiro effect in a superconducting Josephson junction is related to a finite DC 
component of the tunnel current at the resonance frequencies. A similar effect is present 
in a bosonic Josephson junction [Ij. However, current cannot be measured directly and 
furthermore the reservoirs, consisting of the atoms in the left or right well, are finite. 
Therefore, the 'current' changes its sign whenever one reservoir is empty, and in such a 
manner the atoms oscillate between both wells. 

The initial configuration for the following investigations consists of all atoms in the 
lower well. We can characterize the Shapiro resonances by a time average of zero atom 
number difference between the wells. The time averaged imbalance is then, similar as 
in [2], 

(J.)t = ^_^ dt{J.)it). (10) 

A value of 0.5 then corresponds to no population transfer at all, whereas a value close 
to zero indicates a resonance. 

^ We simulated MCTDHB with 4 orbitals [33l [34] and found that 2 orbitals allow for a very good 
description of the physics discussed in this work. 
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3. Shapiro resonances in absence of interactions 

We now discuss the emergence of Shapiro resonances in the atomchip geometry (i.e., by 
driving the double weU separation in the presence of a fixed potential gradient) for the 
non-interacting case. This system can well be captured within the TM model. 

3.1. Enhanced Shapiro effect 

The time-dependence of the control parameter is chosen as 

\{t) = Ao + Ai sin(cjt) , (11) 

with driving frequency uj. In the following we take Aq = 0.675, corresponding to a 
splitting distance of approximately 1 /xm, and a gradient g = 1.5042, corresponding to 
A£'o ^ 271 ' 280 Hz. Due to the linear relationship between A and the tilt AE at Aq 
(see figure [T]), we have to a very good approximation AE{t) = AEq + b • sin(cL;t), where 
b is the driving amplitude of the bias energy. The tunnel coupling instead depends not 
linearly on the control at Aq, but rather in a polynomial fashion. Thus, the general form 
of the tunnel coupling is 

n{t) = no + Y^ n^r^ sin {imcjt) (m = 1, 2, 3, ...) . (12) 



We now decompose the Hamiltonian, equation [2] for Uq = 0, into H = HQ{t) + i^i(t), 
with 

Ho{t) = [nu + 6sin(cjt)] J^, Hi{t) = {AEq - nu))J, - n{t)X • (13) 

Hereby, n is an integer corresponding to the order of the resonance. Next we transform 
into an interaction picture [37] with respect to Ho{t): 

H['\t) = {AEo - nuj)J, - n{t) e-^[""^+^^^^("^)l ^ \k){k\X\k + l){k + 1| 

k 

k 

where we exploited that Jx couples only neighboring states \k) and \k ± 1). Then we 
insert the generating function of the ordinary Bessel functions Ji{z) 



^iz cos (out) 

and obtain 



= j;(.e-*yj,(^), (15) 



Hi'^ = {AEo - nuj)J, - nfj^^ , (16) 
with an effective tunnel coupling ^ 



m=l k=0 
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where Ai = Ji{^' Thus, the contributions to the effective tunnel couphng QP^ add up 

in a vectorial fashion. Further, we have J^^^ = cos (0^) J^ + sin {(/)n)Jy where the phase 0^ 
has no influence on the considered Shapiro dynamics of {Jz){t) when initially all atoms 
are in the lower well. 

We neglected fast oscillating terms and kept only the DC contribution, which is 
reasonable close to the resonance uj ^ A£'o/n, and whenever the driving frequency uj is 
much larger than the coupling Qq [38l [H] . 

Thus, the driven system can be approximately described as an undriven system, 
but with renormalized tunnel coupling [39l[2l|38]. At the resonance condition nu) = AE^ 
there is no energy bias between the left and the right condensate, and thus the atoms 
tunnel at rate ("Shapiro current"). For example, the case n = corresponds to the 
trivial resonance of an unbiased system {AEq = 0) at = 0, with Lorentzian line shape 
characterized by . The eflFect of the driving is then to modify the depth and the 
width of the Lorentzian [39]. A system with flnite energy bias (A£'o ^ 0) has similar 
Lorentzian shaped resonances, but at discrete frequencies uj = AEq/u. 

Most importantly is the oscillatory part of which gives additional cross- 

contributions to QP^ ^ proportional to J^-^^^^ Since > Jn{x) for small values 

of X, those additional contributions to Vt^^ can be considerably larger than those due 
to the flrst term of equation 17, proportional to Jn{^) and these contributions can 
drastically enhance width anaSepth of the resonances. 

A typical resonance structure is reported in flgure [2] (solid line) , showing distinct 
resonances up to order n = 5. We have chosen a driving amplitude which increases with 
the order of the resonance as Ai = 0.03 • AEq/u)^ such that higher resonances become 
more distinct. The n = 1 resonance has not only contributions from -^i^;^^, but also 

from Jo since l^i^^ ^ 0. The n = 2 resonance has not only contributions from 

^2^^)' ^^^^ from -^1(^5) ^i^^^ ^1^^ 7^ from -^o^^) ^^^^^ ^1^^ 7^ 0- Because 

f^^^^ is small for m > 2, the width of the resonances becomes smaller for higher n. 

In the flgure we also compare to the case of an artiflcially constant Q{t) (f]^^^ = 
for m = 1, 2, 3...) (dashed line), which shows only extremely weak resonances for n > 1 
for the same driving amplitude. This is because the effective tunnel coupling is 
proportional to the n-th Bessel function Jni-)^ and our driving amplitude is relatively 



small as typically b ^ n/10. Thus, the contribution of Bessel functions with n > 1 is 
very small. 

For the case of an artiflcially constant AE{t) {b = 0) (dashed-dotted line), we flnd 
very pronounced n = 1 and n = 2 resonances. These resonances are "trivial" in the 



sense that they are due to a cancellation of the oscillatory terms in equation [14} The 
strength of the n-th resonance thus corresponds to the magnitude of which is 

largest for m = 1 and m = 2. 

Hence, the combined driving with AE{t) and Q{t) yields more than the sum of 
driving with only one of them, and we thus term these resonances enhanced Shapiro 
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Figure 2. Shapiro resonances for the non-interacting case. The time- averaged 
population imbalance {Jz)t is plotted, with T = 100. The driving amplitude Ai 
depends non-linearly on the driving frequency cj as Ai = 0.03 • AEq/uj. Solid line: 
Enhanced Shapiro resonances. Dashed line: We compare to the standard Shapiro 
effect (oscillating AE{t), but constant Q). Dashed-dotted line: oscillating but 
constant AE. 



resonances. 



3.2. Optimal choice of parameters 

We shortly discuss the regime of parameters which are best suited to find clear Shapiro 
resonances, taking into account limitations posed by an experiment realization. 



3.2.1. Tunnel coupling Q The tunnel coupling which is determined by the mean 
double well separation, controls how fast the atoms tunnel. Therefore it should be 
large enough in order that the averaging period for equation [TO] is not too long 
(typically QT = 15). If Q is too large, the resonances are shifted according to 
AE VAE'^ + Q'^ [2j. Due to the larger amplitude Rabi oscillations, the resonances 
get wider and less clear. Another disadvantage of too large Q is that for an interacting 
system, the initial state will not be localized [40j . This could be compensated by choosing 
a larger tilt, which has however other disadvantages as discussed below. 



3.2.2. Asymmetry A too small asymmetry (linear gradient g) is unfavorable in 
experimental realizations as it leads to tunnelling dynamics already during the system 
preparation. Furthermore the initial state is not well localized, and no clear resonances 
can be identified when averaging over time. A too large asymmetry on the contrary 
requires large driving amplitudes in order to induce a Josephson current. 
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Figure 3. Shapiro resonances for interactions (a,d) UqN = 1, (b) UqN = 2, and (c) 
UoN = 4 {N = 100). We compare MCTDHB, GP and TM model, (a) We show also 
the case of constant Q or constant AE^ similar as in figure |2j Ai = 0.035 • AEq/uj is 
used, (b) Ai = 0.03 • AEq/uj is used, (c) Ai = 0.025 • AEq/u is used, (d) Scan over 
various Ai. 



3.2.3. Driving amplitude b We have already seen in section [3] that a too small driving 
amplitude reduces the Shapiro current due to the Bessel function structure. A too large 
6, however, has the effect of broadening the resonances Jn{h/(jj)). 

4. Shapiro resonances in presence of interactions 

Resonance shifts 

Shapiro resonances in a many-body system are shifted in frequency due to atom- 
atom interactions. For weak interactions, as in [2], the shift is roughly given by the 
interaction energy [/qA^, with the exact resonance frequencies given in terms of elliptic 
functions [411 E]- However, for typical experimental interaction strengths [/qA > 1 
as we consider here, the resonance frequency depends itself on the driving amplitude 

Ai [391 sg. 

The shift of the resonances due to interactions is demonstrated using MCTDHB 
calculations (solid lines in Figs. [3] (a-c)), where {Jz)t is shown for several interaction 
strengths t/oA = 1, t/oA = 2, and t/oA = 4, respectively. With increasing interaction 
strength, the resonances become shifted towards each other. For relatively strong 
interactions U^N = 4 (figure Is] (c)), the initial state is not completely localized and 
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the contrast of the resonances is reduced. 

The dependence of the resonance shift on the amphtude Ai is shown in figure [3] (d) 
for U{)N = 1, calculated from the TM model. The shift becomes smaller for larger Ai. 
The black line corresponds to Ai = 0.035 • AEq/uj^ which was used in Figs. |3] (a-c) to 
map out the resonance most clearly. 

For UqN = 1 and for the TM model, we compare also to the case of artificially 
driving only with Q{t) or AE{t)^ similar as in figure [2| We find that the enhanced 
Shapiro resonances are much more distinct than the usual Shapiro resonances, and 
the difference is much larger than in the non-interacting case. We also find that the 
resonance shifts of the usual Shapiro resonances (i.e., at constant are different from 
those of the enhanced Shapiro resonances and the case of constant AE. This is similar 
to the findings in [42j, where it has been shown that an oscillating Q{t) leads to a shift 
in the effective interaction strength. 

4-2. Spatial dynamics and damping 

Finally we discuss the differences between several models for calculating the resonances. 
The best description is provided by the MCTDHB method, such that the MCTDHB 
results in Figs. [3] (a-c) (solid lines) serve as a reference. The main question we want to 
answer here are the limitations of simpler models in describing the driven system. 

We compare in Figs. [3] the resonance structure for (a) UqN = 1, (b) UqN = 2, and 
(c) UqN = 4 calculated within different models. For UqN = 1 we find that there is a 
very good agreement between MCTDHB (solid line), GP (dashed-dotted line) and TM 
model (dashed line), a discrepancy only exists regarding the widths of the resonances. 
For UqN = 2 the GP fails to be a good description, especially for the first resonance. 
For UqN = 4: finally, the TM model does neither correctly reproduce depth and location 
of the resonances. The GP equation gives slightly better results than the TM-model, 
but still appreciably misses the correct results. From these simulations we see that, 
whereas for UqN = 1 and UqN = 2 the coincidence between TM and MCTDHB is very 
good for all times, for UqN = 4 the coincidence is lost at short times. Thus, whereas for 
weaker interactions Q{t) and AE{t) as deduced from static orbit als are sufficient, the 
full self-consistency of MCTDHB is needed for stronger interactions {UqN = 4). 

The spatial dynamics of the condensates as shown in figure [i] (a) is adiabatic 
to some extent, since the driving frequency is much smaller than the transverse trap 
frequency AEq <C (jOx.y However, as we have seen, for strong interactions {UqN = 4) 
the exact details of the condensate oscillations become important and have to be taken 
into account. 

We show in Figs.[4](b-c) example trajectories for the first resonances of (a) UqN = 1, 
(b) UqN = 2, and (c) UqN = 4, with same color code as in figure |3| The upper panels 
show the atom number imbalance (J^)(t), whereas the middle panels show the time 
averaged atom number imbalance {Jz)t = f Sq ^t' {Jz{t')) . The trajectories of the TM 
model and MCTDHB mostly deviate, in particular for UqN = 4, where the trajectory 
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Figure 4. (a) Density for UqN = 1 and driving frequency u/AEq = 0.85. {N = 100) 
(b-d) Example trajectories for different interaction strengths {N = 100), with color 
code as in figure [sj The upper panel represents {Jz)(t)^ where the gray band represents 
the variance /S.Jz for MCTDHB, the middle panel the averaged quantity {Jz)ti and 
the lower panel the fragmentation (pn — p22)/N. (b) UqN = 1 and uj/AEq = 0.85. (c) 
UqN = 2 and uj/AEq = 0.6. (c) UqN = 4 and uj/AEq = 0.52. The gray line represents 
an MCTDHB solution for N = 1000 (with UqN constant). 



of the TM model differs completely from the MCTDHB one after 10 ms. In order to 
point out the importance of the renormalization to the tunnel coupling according to the 
last term in equation [Tj we also indicate results where it is not included (bright dashed 
line figure [3] (c)). 

The GP equation is in principle not a good description of the Shapiro dynamics, 
as it does not take into account damping of the Shapiro oscillations on a time scale 
of tens of ms. We show in the upper panels of Figs. [4] (b-c) the width of the number 
distribution AJ^(t) for the TM calculations by the bright band around the mean value 
{Jz){t). Whereas for short times the number distribution is very narrow, it starts to 
become much broader after a time which depends on interactions. At the same time 
the GP results start to deviate from the exact ones. Finally, the distribution stays 
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very broad and the mean value does not really change anymore, while the GP equation 
predicts oscillations of {Jz){t) with approximately constant amplitude. 

Thus, the time averaged imbalance {Jz)t does not fully characterize the system 
dynamics, and it is more appropriate to study the time evolution of the imbalance. 
However, it has been expressed in the literature [2] that for weak interactions and the 
"standard" Shapiro effect, the averaged imbalance {Jz)t of the TM model and the GP 
equation should be the same, although the dynamics are very different in each case. 
For the configuration we discuss here we find this to be justified only for UqN = 1 (see 
figure |3] (a)). In this case the damping is slower than the period of the oscillations, and 
{Jz)t corresponds approximately to the average of the extremal values, see the middle 
panel of figure [3] (b). However, for stronger interactions, the imbalance becomes stuck 
at an early point of the dynamics and thus the results for the averaged imbalance are 
significantly different from those predicted by the GP equation. 

In order to have another perspective on the dynamics of damping, we plot the 
population difference of the natural orbitals (/>i(x) (i = 1, 2). As can be seen in the lower 
panels of figure [4| during the Shapiro dynamics the BECs become fragmented into two 
incoherent {{d\a2) = 0) condensates due to the nonlinear interactions. We find that 
indeed the GP description starts to fail whenever the system starts to fragment. 

In addition we plot in figure [i] (d) also results for = 1000 atoms (thin line with 
symbols), in order to demonstrate that fragmentation is clearly important also for larger 
numbers of atoms. Thus, convergence of MCTDHB results to GP results happens only 
at very large N. 

Conclusion 

In conclusion, we have studied Shapiro resonances in a configuration where not only the 
bias potential, but also the tunnel coupling is driven dynamically. This is typical for 
double well potentials realized on atomchips and thus our findings are directly relevant 
for future experiments. We show that this configurations has favorable properties as 
it leads to enhanced Shapiro resonances. Due to a spatial deformation of the potential 
induced by the driving, the question of transverse excitations is of great interest. We 
find that, at least for significant interactions, the realistic MCTDHB method has to 
be used instead of simpler models in order to properly capture spatial dynamics of the 
involved modes. 
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